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Abstract 

We consider a model of a square-wave bursting neuron residing in the regime of tonic spiking. Upon 
introduction of small stochastic forcing, the model generates irregular bursting. The statistical properties 
of the emergent bursting patterns are studied in the present work. In particular, we identify two principal 
statistical regimes associated with the noise-induced bursting. In the first case, (type I) bursting oscilla- 
tions are created mainly due to the fluctuations in the fast subsystem. In the alternative scenario, type II 
bursting, the random perturbations in the slow dynamics play a dominant role. We propose two classes 
of randomly perturbed slow-fast systems that realize type I and type II scenarios. For these models, 
we derive the Poincare maps. The analysis of the linearized Poincare maps of the randomly perturbed 
systems explains the distributions of the number of spikes within one burst and reveals their dependence 
on the small and control parameters present in the models. The mathematical analysis of the model 
problems is complemented by the numerical experiments with a generic Hodgkin-Huxley type model of 
a bursting neuron. 



1 Introduction 



Differential equation models of excitable cells often include small random terms to reflect the unresolved or 
poorly understood aspects of the problem or to account for intrinsically stochastic factors lfTll8ll9l [T0l[T5l[T6l 
[32l|4T]|39j|43l|46l. In addition, many neuronal models also exhibit multistability 113811261 . In systems with 
multiple stable states, noise may induce transitions between different attractors in the system dynamics, thus, 
creating qualitatively new dynamical regimes, that are not present in the deterministic system. In the present 
paper, we study this situation for a class of square-wave bursting models of excitable cell membranes. This 
class includes many conductance-based models of excitable cell membranes. Here we just mention the 
model of a pancreatic (3— cell |0[7]], models of neurons in various central pattern generators such as those 
involved in insect locomotion l20ll . control of the heartbeat in a leech |[25l . and respiration in mammals 
El [5), to name a few. These models, as well as the underlying biological systems, exhibit characteristic 
bursting patterns of the voltage time series: clusters of fast spikes alternating with pronounced periods 
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Figure 1: The dynamical patterns generated by a model of of a square- wave bursting neuron (11.11) and (11.21) : 
(a) periodic bursting and (b) tonic spiking. 

of quiescence (Fig. [T^). For introduction to bursting, examples and bibliography, we refer the reader to 
|[26l I3T1 l37l l38l . The dynamical patterns generated by the conductance-based models typically depend 
sensitively on parameters. For example, models of square-wave bursting neurons often exhibit both bursting 
and spiking behaviors for different values of parameters (see Fig. [U,b). In many relevant experiments, 
the transition from spiking to bursting is achieved by changing the injected current. In the present paper, 
we consider a model of a square-wave bursting neuron in the regime of tonic spiking (Fig. \Vp). We show 
that a small noise can transform spiking patterns into irregular (noise-induced) bursting patterns and describe 
two distinct mechanisms for generating noise-induced bursting. In the first scenario, bursting oscillations are 
triggered by the fluctuations in the fast subsystem. We refer to this mechanism as type I bursting. In contrast, 
the bursting dynamics in type II scenario are driven by the random motion along the slow manifold. For 
each of these cases, we describe the statistical properties of the emergent bursting patterns and characterize 
them in terms of the small and control parameters present in the model. 

Noise-induced phenomena have received considerable attention in the context of neuronal modeling (see, 
e.g., flU [8] [32l [39l |4l] |43] @6l). A representative example is given by a 2D excitable system perturbed 
by the white noise of small intensity d. In the presence of noise and under certain general conditions, 
a typical trajectory occasionally leaves the basin of attraction (BA) of the stable equilibrium and makes a 
large excursion in the phase plane of the deterministic system before returning to a small neighborhood 
of the stable fixed point (Fig. |2k). This gives rise to irregular spiking (Fig. |2j)). The properties of the 
noise-induced spiking and stochastic resonance type effects arising in the context of the perturbed FitzHugh- 
Nagumo model have been considered in [fl] [U [U [TOj (see also |f3l [17] (TU [191 for the mathematical analysis 
of more general classes of related phenomena in randomly perturbed slow-fast systems). In the present 
paper, we study a related mechanism for irregular bursting. Specifically, we consider a class of models of 
square-wave bursting neurons: 



where / and g are smooth functions and < e <C 1 is a small parameter. We refer to (11.11 ). where y is 



x = f(x,y), 

y = eg(x, y), x = {x\x 2 ) T G M 2 , y G R 
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Figure 2: (a) A phase-plane trajectory of the randomly perturbed FitzHugh-Nagumo model in excitable 
regime (see HI for the model description and the parameter values), (b) The time series corresponding to 
the phase plot in (a). 

treated as a parameter, as a fast subsystem. It is formally obtained from dl.ll ) and dl.2b by setting e = 0. 
We assume that the fast subsystem has a family of stable limit cycles and that of stable equilibria for y in 
a certain interval y £ (y S n,ybp) (see Fig. [3^)- The additional assumptions on dl.ll ) and (11.2I ). which are 
explained in Section 2, imply that for small e > 0, dl.ll) and (1 1 .2b has a stable limit cycle as shown in Fig. 
|3j;. In the presence of noise, a typical trajectory of the randomly perturbed system will occasionally leave 
the BA of the limit cycle of the deterministic system to make an excursion along the curve of equilibria of 
the degenerate system, E (see Fig. Hk). Thus, in analogy to the 2D FitzHugh-Nagumo model (Fig. [2^), 
noise transforms spiking dynamics into irregular bursting. We refer to the latter as noise-induced bursting. 
In both examples above, irregular spiking (Fig. |2^) or bursting patterns (Fig. @^,b) are created due to the 
escape of a trajectory of the randomly perturbed system from the BA of a stable fixed point in the case 
of spiking or of that of the stable limit cycle in the case of bursting. The statistics of the first exit times 
can then be related to the properties of the emergent firing patterns such as the frequency of spiking or the 
distribution of the number of spikes within one burst. Compared to the analysis of the irregular spiking in 
the randomly perturbed FitzHugh-Nagumo model (Fig. [2b. , the analysis of the noise-induced bursting faces 
several additional challenges due to the fact that in the latter case one has to consider the exit problem for 
the trajectories near a stable limit cycle as opposed to those near a stable equilibrium in the former case. The 
structure of the BA of the limit cycle combined with the slow-fast character of the vector field determines the 
main features of the resultant bursting patterns. The description of the principal statistical regimes associated 
with the noise-induced bursting is the focus of the present paper. 

There are general mathematical approaches for analyzing exit problems for stochastic processes gener- 
ated by randomly perturbed differential equations such as dl.lb and dl.2b : the Wentzell-Freidlin theory of 
large deviations lfl9ll and the geometric theory for randomly perturbed slow-fast systems due to Berglund 
and Gentz [3 ]. In this paper we study the vector fields arising in the context of bursting. The specialized 
structure of this class of problems allows us to keep the analysis of the present paper self-contained and 
avoid using more technical methods, which are necessary for analyzing more general situations. Our analyt- 
ical approach is based on the reduction of a randomly perturbed differential equation model to the Poincare 
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map and studying the exit problems for the trajectories of the discrete system. Using maps is quite natural 
in the context of bursting due to the intrinsic discreteness of bursting patterns imposed by the presence of 
spikes. Reductions to maps have been very useful for analyzing bursting dynamics in a variety of determin- 
istic models (6l[33][34l[35l|40]]. As follows from the results of the present paper, the first return maps also 
provide a very convenient and visual representation for the mechanism underlying noise-induced bursting. 
In particular, we show that the distributions of spikes in one burst in many cases are effectively determined 
by ID linear randomly perturbed maps. We develop a set of probabilistic techniques for analyzing the dy- 
namics of randomly perturbed ID and 2D linear maps such as those arising in the analysis of bursting. The 
special structure of this class of problems, which is motivated by the applications to bursting affords a more 
direct and simpler analysis than the treatment of more general classes of random linear maps found in the 
literature ll29ll2Tll28ll431. 

The outline of the paper is as follows. In section 2, we formulate our assumptions on the deterministic 
system. We then present the preliminary numerical results, motivating our formulation of the randomly per- 
turbed models at the end of this section. Specifically, we distinguish two types of the noise-induced bursting. 
Type I bursting is generated due to the fluctuations predominantly in the fast subsystem, while type II burst- 
ing is induced by variability mainly in the slow variable. Accordingly, we introduce two types of models that 
generate type I and type II bursting patterns. Section 3 develops a set of probabilistic techniques, which will 
be needed for the analysis of the first return maps for the randomly perturbed differential equation models. 
We first analyze a simple linear map with an attracting slope and small additive Gaussian perturbations in 
Section 3.2. Due to the simple structure of the map, we obtain very explicit characterization of the first exit 
times for this problem. The analysis of this first relatively simple example provides the guidelines for the 
more complex cases dealt in Sections 3.3-3.5. Section 4 contains the definition and the construction of the 
Poincare map for the type I randomly perturbed model introduced in Section 2. The 2D Poincare map is 
decomposed into two ID maps for the fast and slow subsystems, which are constructed in Sections 4.2 and 
4.3 respectively. In Section 4.4, we apply the results of Section 3 to the linearization of the Poincare map to 
derive the distributions of the first exit times. The latter are interpreted as the distributions of the number of 
spikes in one burst. In Sections 4.5, we outline the modifications necessary to cover type II models. Since 
the analysis for type II models closely follows the lines of that for type I models, we omit most of the details. 
Finally, the numerical experiments in Section 5 are designed to illustrate our theory. 



2 The model 

In the present section, we introduce the model to be studied in the remainder of this paper. We start by 
formulating our assumptions on the deterministic model and then describe the random perturbation. 

2.1 The deterministic model 

We consider slow-fast system (11.11 ) and (11.21 ) in M 3 with one slow variable. The fast subsystem associated 
with (11.11 ) and (11.21 ) is obtained by sending e — > in (11.21 ) and treating y as a parameter: 

x = f(x,y). (2.1) 
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Figure 3: (a) The bifurcation diagram of the fast subsystem (12.11) . L denotes a cylinder foliated by the 
stable periodic orbits. The lower branch of the parabolic curve E is composed of stable equilibria of the fast 
subsystem (see Fig. [6}) for the plot of a representative phase plane of the fast subsystem for y G (y S m ybp))- 
(b,c) Periodic trajectories of the full system (11.11 ) and (11.21 ) are superimposed on the bifurcation diagram of 
the fast subsystem. Assumptions (SE) and (SB) (see the text) result in a bursting limit cycle plotted in red 
in (b), while (SS) yields spiking (c). 

Under the variation of y, the fast subsystem has the bifurcation structure as shown schematically in Fig. 
[3^. Specifically, we rely on the following assumptions: 

(PO) There exists y^p G R such that for each y < yt p , Equation (12.11 ) has an exponentially stable limit cycle 
of period T(y): 

L(y) = {x = <P(s, y): < s < T(y)}. (2.2) 
The family of the limit cycles, L = \J y<yb L(y), forms a cylinder in M 3 (Fig. |3^). 

(EQ) There is a branch of asymptotically stable equilibria of d2.ll ), E = {x = ip(y) : y > y S n}, which 
terminates at a saddle-node bifurcation at y = y sn < yt p (Figured). 

(LS) For each y G R, the cu— limit set of almost all trajectories of d2.il ) belongs to L(y) {J{ip{y)}- 

Remark 2.1. At y = y^, L, either terminates or + 0) looses stability. We do not specify the type 

of the bifurcation at y = yb p - It may be, for instance, a homoclinic bifurcation as shown in Fig. [3^, or a 
saddle-node bifurcation of limit cycles 1221 . 

Having specified the assumptions on the bifurcation structure of the fast subsystem, we turn to the slow 
dynamics. The geometric theory for singularly perturbed systems implies the existence of the exponentially 
stable locally invariant manifolds E e and L e , which are O(e) close to Ef]{(x,y) : y > y sn + 6} and 
Lf]{(x,y) : y < yt, p — 5}, respectively, for arbitrary fixed 5 > and sufficiently small e > lfl4l l27l . 
Manifolds E e and L e are called slow manifolds. For small e > 0, the dynamics of dl- lb and (11.21 ) on the 
slow manifolds is approximated by 

L e - y = eG{y), y < y bp - 5, (2.3) 

Ee- y = eg(ip(y), y), y>y sn + 5, (2.4) 
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where 



1 rT(y) 



(2.5) 



We distinguish two types of the asymptotic behavior of solutions of (11.11 ) and (11.21 ): bursting and spiking 
(see Fig. [T|). The following conditions on the slow subsystem yield bursting. 

For some c > independent of e, 
(SE) 

9&(y),y) < ~c for y>y sn , (2.6) 

(SB) 

G(y) > c for y < y bp . (2.7) 

Under these assumptions, for sufficiently small e > a typical trajectory of (11.11 ) and (11.21 ) consists of 
the alternating segments closely following L e and E e and fast transitions between them (see Fig. [3}}). 
For detailed discussions of the geometric construction of 'bursting' periodic orbits, we refer the reader to 
II3T1I3711 . To obtain spiking, we substitute (SB) with 

(SS) G(y) has a unique simple zero at y = y c G (y sn , ybp): 

G(y c ) = and G'(y c ) < 0. (2.8) 

In this case, the asymptotic behavior of solutions follows from the following theorem due to Pontryagin and 
Rodygin: 

Theorem 2.2. H36\l lfe>0is sufficiently small, di.il ) and rti.2D has a unique exponentially stable limit cycle 
L e (y c ) of period T{y c ) + 0(e) lying in an 0(e) neighborhood of L(y c ), provided (SS) holds. 

Almost all trajectories of (11.11 ) and (11.21 ) are attracted by the limit cycle lying in an O(e) neighborhood 
of L(y c ). This mode of behavior is called spiking (see Fig. [3fc and Fig. [3})). In the remainder of this paper 
we assume (SS), in addition, to (PO), (EQ), (LS), and (SE). 

2.2 The randomly perturbed models 

In this subsection, we provide a heuristic description of the effects of the random perturbations on the 
dynamics of (11.11 ) and (1 1 -2b - To study these effects quantitatively, at the end of this section, we propose two 
randomly perturbed models. 

Suppose the trajectories of (11.11 ) and (11.21 ) experience weak stochastic forcing, such that the perturbed 
trajectories represent well-defined stochastic processes and are close to the trajectories of (11.11 ) and (11.21) on 
finite intervals of time. Since the trajectories of the unperturbed system remain in a small neighborhood 
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Figure 4: Noise-induced bursting, (a) A trajectory of the randomly perturbed system is shown in the phase 
space of the frozen system (11.11 ). (11.21 ) with e = 0. The trajectory leaves the basin of L(y c ) mainly due to 
the fluctuations in the fast plane. This is characteristic to type I bursting. An alternative type II scenario is 
shown in plot (b), where the fluctuations in the slow direction dominate in the mechanism of escape from the 
basin of the stable limit cycle. The trajectory in (b) samples a wide region of L and leaves a neighborhood 
of L near the right boundary, y « while that in (a) remains near L(y c ) most of the time and jumps down 
near y as y c . The differences translate into the distinctive features of the generic time series of the bursting 
patterns generated via type I or type II mechanisms shown in plots (c) and (d) respectively. Note that the 
longer burst in (c) has a typical square-wave form (roughly, determined by L{y c )), while the burst shown in 
(d) exhibits more variability due to the drifting of the trajectory along L. 
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of L{y c ) (possibly after short transients), we expect that in the presence of noise the trajectories will occa- 
sionally leave the BA of L(y c ) and after making a brief excursion along E will return back to the vicinity 
of L(y c ). Therefore, under random perturbation the system can exhibit bursting dynamics, while the un- 
derlying deterministic system is in the spiking regime. We refer to this mode of behavior as noise-induced 
bursting. Our goal is to describe typical statistical regimes associated with the noise-induced bursting and to 
relate them to the structure of (11.11 ) and (11.21) and to the properties of the stochastic forcing. To illustrate the 
implications of the structure of the deterministic vector field for the bursting patterns that it produces under 
random perturbations, we refer to the following numerical examples. Note that the BA of L(y c ) naturally 
extends along the cylinder of periodic orbits L (Fig. [3fc). The escape from the BA of L(y c ) can be dominated 
by the fluctuations along L or by those in the transverse plane. These two possibilities are shown in Fig. [4] 
The trajectory shown in Fig. @k spends most of the time near L(y c ) and leaves its BA due to the fluctuations 
in the fast subsystem. We refer to this scenario as type I escape. Alternatively, the trajectory shown in Fig. 
HJ) travels a good deal along L before the escape and exits from the BA near y = y^. This mechanism is 
dominated by the slow dynamics. We refer to this scenario as type II escape. These mechanisms of escape 
translate into distinct features of the resultant bursting patterns. First, note that since in type I and type II 
scenarios, the transition from spiking to quiescence typically takes place at y y c and y « yb p respectively, 
by (11.21) and (EQ), the corresponding interburst intervals are approximately equal to 



In addition, we expect that the interspike intervals (ISIs) within one burst in type I scenario are localized 
about T(y c ), since the trajectory of the randomly perturbed system in the active phase of bursting spends 
most of the time near L(y c ). In type II bursting patterns, ISIs are expected to have more variability, since 
the trajectories sample a wider range of ISIs during their excursions along L. Perhaps, a more pronounced 
distinction between these two types of bursting patterns lies in the degree of the variability of the spikes 
in one burst. Most of the spikes forming a burst in type I pattern are generated by (12.11 ) with y as y c and, 
therefore, are similar in shape (Fig. @fc). In contrast, spikes in type II scenario are subject to more variability 
and the bursting patterns typically have ragged shape (Fig. [4}i). 

To study type I and type II noise-induced bursting patterns it is convenient to consider two types of 
models. Type I model incorporates random forcing in the fast subsystem: 




where 




type I, 
type II. 



x t = f (x t ,y t ) + apw t , 
yt = eg(x t ,y t ), 



(2-9) 
(2.10) 



while, in type II model the slow subsystem is forced 



xt = f{x t ,y t ), (2.11) 
y t = e(g(x t ,y t ) +o-qw t ) . (2.12) 

Here, < a -C 1, p(x, y) = (p 1 (x, y),p 2 (x, y)) T and q(x, y) are differentiable functions; wt stands for the 



white noise, i.e. a generalized derivative of the Wiener process. 
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3 The randomly perturbed maps 



In this section, we develop probabilistic tools needed for the analysis of randomly perturbed systems (12.91 )- 
(12- 12b - The number of spikes in one burst is a natural random variable associated with the noise-induced 
bursting. It is commonly used in the experimental studies of bursting and we shall adopt it for characterizing 
irregular bursting patterns in this work. In Section 4, we will show that the number of spikes in one burst is 
represented by a stopping time (more precisely, the level exceedance time) of a discrete random process, the 
Poincare map of the randomly perturbed system (|2.9I >- (I2.12I ). In preparation for the analysis of the linearized 
Poincare map in Section 4, in the present section we study certain stochastic linear difference equations. The 
equations of this form equations have been considered in the literature before. The study was initiated by 
Kesten |[29l who considered multidimensional case (in which the coefficients of the stochastic equations 
are random matrices). Subsequent work focused mostly on the ID case. We refer the reader to the papers 
|2T1 l45l . which contain representative results, examples of applications, and further references. There is 
also a review paper lfl2l . unfortunately not easily accessible. The convergence properties of the solutions 
that we will need could be deduced from a general theory of stochastic difference equations. However, the 
results in the literature are often stated in the most general form and some of the proofs are rather involved. 
We will be dealing with special cases that are much easier to justify. For this reason, and also to keep the 
paper self-contained we will include the proofs of the needed results. 



3.1 Geometric random variables 



We begin by recalling the necessary properties of geometric random variables (RVs). Recall that Y is a 
geometric RV with parameter p, < p < 1 if 

F(Y = k) =p(l -pf" 1 , k>l. (3.1) 

We refer the reader to ll28l Chapter 5] for the review of the properties of geometric distributions and their 
applications. In particular, the following characterization of geometric RVs is classical. 

Lemma 3.1. Let Y be a RV with values in the set of positive integers. Y is a geometric with parameter p, 
0<p<l,iff 

P(y = n) = pF(Y >n), n > 1. (3.2) 



Lemma [37TI motivates the following definition: 

Definition 3.2. Let Y be a random variable with values in the set of positive integers and let < p < 1. We 
say that Y is asymptotically geometric with parameter p if 

F(Y = n) 

lim f = p. (3.3) 

n^oo F(Y > n ) F v 
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3.2 The randomly perturbed map: additive perturbation 



Consider 

Y n = Ay n _i+?r n , n>l, (3.4) 

where r\ , r2 , ■ ■ ■ are independent identically distributed (IID) copies of the standard normal RV, and Yq is 
a real number. We will use N(/j,, rj 2 ) notation for a normal RV with mean /i, variance rj 2 , and probability 
density function given by 

i r (x-v) 2 \ 

, exp < ^ — > , — oo < x < oo. 

We will also let Z denote a generic N(0, 1) RV and we will write 

X 



$( S ) := -L /" e- i2 / 2 dt, 

y 



for its distribution function. For a given h > 0, let 

t = inf{A; > 1 : Y k > h}. 

Theorem 3.3. Let 

EG (0,1), A = l-e, /3 2 = — ^ r, and /i-y >0. (3.5) 

e(2 — ej 

Then for sufficiently small q > 0, r « asymptotically geometric RV with parameter 

We precede the proof of the theorem with the auxiliary 
Lemma 3.4. For n > l,Y n is a normal RV with 

EY„ = X n Y and var Y n = — \ =: 0i. (3.7) 

1 — X A 



In particular, 



Y n Y = N(0,l3 2 ), 



where — —* (and =) denote the convergence (equality) in distribution. 

Proof (Lemma [3~4T >: The statements in (13.71 ) are verified by a straightforward calculation. The rest follows, 
because EF n ^0 and (3 n — * (3. 

Proof (Theorem [331): Let Y fc * = max{Y, : 1 < j < k}, k > 1. Then 

p( T = n + 1) = P(y n+1 > y: <h) = p(y n+1 > h|Y n * < h)¥(Y* < h) 

= P(y n+ i > h\Y n < h, y n _! <h,...,Y < h)¥(r > n + 1) 

= P(Y n+1 >/i|Y n <fc)P(r>n + l). (3.8) 
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In the last equality, we used the fact that {Y n } is a Markov process which is clear from (13 .4b - By (13 - 8b . 

p ._ P(r = n + l) _ P(y B+1 >^,y B <fe) 
Pn -"P(r>n + l) " riyn+1>A|yn - ftJ " P(y„</i) • ( j 

In accordance with Definition 13.21 we need to show that {p n } converges and to estimate the limit. By 
Lemma [3~4l 

P (Y n <h) — > &(h/P), as n -» oo. 
Next, we turn to estimating the numerator in ( I3.9I ). We have 

Q n := P (y n+ i > h, Y n < h) = P (XY n + qr n+l > h,Y n < h) 
-» P(AY + > ft,y < /i) =: Q, 

where Z is standard normal, Y is iV(0,/3 2 ) and they are independent. This follows from Lemma [3~4l and 
the fact that r n+ \ is N(0, 1) and is independent of Y n . Q is the probability that a 2D Gaussian vector is in 
the region [h, oo) x (— oo, h]. There are several ways of estimating this probability. We take the following, 
elementary approach. Let X = h — Y so that X is iV (h, /3 2 ) and is independent of Z. Then 

e, l-e„ „ A 1 f^^f^ eh+(l-e)s\ =is=*£ 



Q = ¥[Z>-h + - — -X,X>0}=^^l P ( Z > ' v ^ — ^le^^ds. 



By the well-known asymptotics (see 11131 Ch. VII, Lemma 2 and Sec. 7, Problem 1]) 



u 2 

1 e~ir / „ / 1 



P(Z > «) = 1 -$(«) = —=——( 1 + ( — )) , n>0. (3.10) 



Hence, for sufficiently small <j > <C e), we have 



i 1 Ae/i+(l-e)s) 2 . (s-h) 2 \\ 



Since 



we obtain 



Q-iTTj " -TTTi \ (3.11) 

27T (5 J eh + (l-e)s 



(eh + (l-e)s) 2 (s-h) 2 _ {s-ehf h 2 



c 2 p 2 c 2 /? 2 ' 



2vr/3 6XP j 2/3 2 J 7 e/i + (1 - e) s ' 



By Laplace's method B71 . for sufficiently small q > (<j <C e), the last integral is asymptotic to 

V2n V2tt<; 



(hs+(l-s)sh)^T/^ he(2-eY 

Hence, 



2^/^(2 -e) eXP i 2/52 / 2/3 2 
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By the same reasoning the error term from (13.10b is of order 

h 2 } 



which gives ( 13.61 ). □ 



3.3 The randomly perturbed map: random slope 

Consider a process 

Y n = p(l + 0Ti in )Y n _i + 0T 2 ,n, n > 1, (3.12) 

where (ri n , T 2jn )^L 1 are IID copies of a two dimensional random vector (r!,r 2 ). Here, we assume that 
(ri,r 2 ) has bivariate normal distribution with mean vector and covariance matrix S 2 = [ojj], where 
= cov(rj, rj), 1 < i, j < 2. We assume that the entries <Tj j are of order 1 in a sense that they do not 
depend on other parameters. Recall that the probability density function of a multivariate normal random 
vector (ri, . . . , r^) with mean vector and covariance matrix £ is given by 

exp < — x X - x > , x = (xi, . . . , Xrf) . 



• v /(27r) d det(E) I 2 
and we denote such vectors by iV(0, E). 
For a given h > 0, let 

t = inf{£; > 1 : 1*. > h}. 

Theorem 3.5. Suppose that h and /i £ (0, 1) are both of order 1 and a <C 1 so that the following condition 
holds 

7 := pE\l + ari\ < 1. (3.13) 
Then r is asymptotically geometric RV with parameter 

p=^=e-£(l + 0(a 2 )), (3.14) 



where a positive constant c depends on h, p, and S 2 , but not on a. 

As before, we first establish convergence of {Y n } and characterize the limit. Iteration of (13.121 ) yields 

Y n = fi(l + OTi in )y n _i + ar 2 , n = Ml + <?ri,n) (M 1 + <m,Ti-l)^-2 + ^2,n-l) + <7^2,n 
n n— 1 n 

= •••=^ n yon( i + fjr ^)+ f7 Z]^ r ' 2 .«-i n c i + <7r i^), (3.15) 

3=1 3=0 <=n-3'+l 

where as usually, JI j=fc (*) = lifA;>m. 
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Lemma 3.6. 



3-1 



Y n — > Y = a ^2 ^92,j + (*9x,i), n -> oo, 

j=0 £=0 



(3.16) 



where (gij,g2,j), j = 0,1,2,... are 7/D copies of two-dimensional random vector, which is equal in 
distribution to (ri,r 2 ). 



Proof (Lemma 13 . 6l > : First, we show that Y is well-defined as the series in (13.161 ) converges almost surely. To 
see this, note that the summands 

3~1 

g2,i\\{i + <?gi,t) 

1=0 

are martingale differences with respect to the natural filtration. By triangle inequality, independence, and 
(13131 . 



E 



m j—1 
j=0 1=0 



<aE\g 2 \J2^ 



3=0 



3-1 



e=o 



aE\g 2 \ jr J (E|l + ortf = ^(1 - 7 (m+1) ) < 

3=0 ' ' 



Hence, the partial sums of the right-hand side of (13.161) form an Li-bounded martingale which converges 
almost surely by the martingale convergence theorem (see e.g. ll42l ). For every n > 1 



n-l 



n-1 



3-1 



i=0 €=n-j"+l i=0 £=0 

Since the sequence on the right converges almost surely and the almost sure convergence implies con- 
vergence in distribution, we infer that the sequence on the left converges in distribution. To conclude that 

Y n — > Y it is enough to show that the first term on the right-hand side of ( 13.151) converges to in probability. 
But that is clear since we have 



E 



Y of i n l[(l + ar K 

3=1 



|y |/x™n iE l 1 + fTri jl = l y ol^- 

3=1 



Hence, by Markov inequality it goes to in probability. 



□ 



Proof (Theorem [331 ): The proof follows the lines of the proof of Theorem 13.31 The main complication in 
treating the present case is that we know less about the distribution of Y n than in before. Nonetheless, we 
will argue that for large n 



Pn ■-- 



ir(r = n) 
P(r > n) 



P(//(l + <rri in )y n _i + o-r 2 , n > h\Y n -i < h) 



(3.17) 
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is approximately constant. For this, we rewrite the right hand side of (13.171 ) as 

PQu(l + ar^Yn-i + ar 2i n > Mn-i < h) 

and since the denominator converges to F(Y < h) we focus on the numerator. Let (ri,r 2 ) be a generic 
vector distributed like (ri >n , r 2jn ) and independent of Y . Since for every n > 1, (n,nj r 2,n) is independent 
of y n _i, as m oo we have 

(n,fi,^2,n>^n-l) — >• (n,r 2 ,y). 

Thus, 

P(/x(l + ari)y„_i + ar 2 > h, Y n ^ < h) — ► P(/i(l + ar x )Y + ar 2 > h,Y < h), 
which establishes the existence of p = lim n ^oo p n . 

To estimate p, we first recall that (n, r 2 ) is bivariate normal if and only if every linear combination of 
r\ and r 2 is a normal RV. Hence, conditionally on Y = y, a(fiyri + r 2 ) is iV(0, o" 2 o"y) RV, where 

oj = 4j + mV^ii + 2 W (Ji 2 . (3.18) 

Therefore, 

P(/x(l + <7ri)y + <rr 2 > /i, Y" < h) = F(a{fiYn |r 2 )>li-^,l'< h) 

= I" HZ > h -^)dF Y{ y) = I" (l - * ( "-M)) dF Y (y) 
J-oo VOy J-oo\ V Ga y J J 

l-$( h -M°))F(Y<h), 



where — oo < y < h by the mean value theorem. Hence, 

= P(/x(l + <"-i)y + OT- 2 > fa, y < h) = ( h-nyp 
P P(Y<h) V ™vo 

Let c := c(yo) where 

h — fix h — fix 



c(x) = c h ^T, 2 (x) :-- 



o x y 7 ^v'tiX 2 + IfiOnX + cr 2 2 

Then, by (13101) 



c\ cr „ / cr 2 



p = 1 - * - = — = e"^ 1 + 



Furthermore, by elementary analysis we see that: 



c(x) is increasing on x E (— oo, x*) and decreasing on x G (x*, cxd), where 

* CJ 2 ! + /l(Ti 2 



fi{ho$ 2 + 0-12) 
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(l-M)fc 



C( OO) CT 11 , C{h) (( Mh£ril )2 + 2^ ri2 / l+0 . a2 )V2 (( M / lCTll +a 22 )2-2 At / l ( ( 7 11 a 22 -a 12 )) 1 / 

given by a quite unwieldy expression that depends on /t and S2 but not on 



T72, and c(x*) is 



In particular, c is bounded away from and 00 provided p, and h are positive and \i < 1. This proves (13.141 ). 

□ 



3.4 A two-dimensional randomly perturbed map 

In this subsection we consider the following two dimensional model: 

Cn+l = fJ-Cn (1 + OTi,rH-l) + <™"2,n+l> 

?7n+i = A?7„ + ecrr3 „+i + ea2^«- 



(3.19) 
(3.20) 



where (ri, n , r2,n, r^ n ), n > 1, is a sequence of IID copies of (ri,r2,rs) which, as follows form adiscussion 
at the beginning of Section l4~4l is assumed to be a trivariate normal random vector N(0, S3), with S3 = 
1 < i>i 3 < 3, where o-jj = cov(ri,rj) do not depend on any parameters in (14.441 ) and (14.45b - For 
positive h\,}i2 = 0(1), we define 

r ? = inf{£ fc > t„ = inf{r/ fc > /i2}- 

fc>l fc>i 

We are interested in r = min{r^, t„}. We know the distribution of rg from Theorem l3.5l As we will show 
below, under the suitable conditions the distribution of r is again asymptotically geometric. Moreover, if 
e > is small then t v has practically no effect on the distribution of r. 



In order to be more precise, let us define 

_ /i(l + OTl )n ) 

ea 2 A 



1 G n — 


r2,n 







and G n 



Then, (l3~T9l) and (13201) are described by 

8^+1 = A n+1 6 n + o-Gn+i, n>l. 



(3.21) 



(3.22) 



Theorem 3.7. Le? /u, a, e E (0, 1) fte smc/i ?/ia? /x of order 1 and a <C 1 *o condition H3.13i holds. 
Assume e -C larccf sef A = 1 — e. Suppose further that h\ and h-2 are of order I. Then r is approximately 
geometric RV with parameter p satisfying 



a 



cV2vr 



e 2d- 3 



(3.23) 



a«(i where the constant c depends on h\, p, and S3 but not on a. 

The following lemma shows that {G n } converges in distribution and describes the limit. 
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Lemma 3.8. 

Q n ^X±aY, ) G k> n^oo, (3.24) 

k=i \j=l J 

where A n and G n ,n = 1,2, dots are defined in A3. 21 I ). Furthermore, this random vector X satisfies the 
distributional equation 



X = AX + o-G, 



where 



A 



li(l + an) 
ea 2 A 



and G 



r-2 
er 3 



(3.25) 



(3.26) 



(ri, r 2 , T"3) iV(0, £3) fte generic copies of A n and G n , and, X on the right hand side of A3. 251 ) is indepen- 
dent of (A, G). 

Proof (Lemma I3T8T ) : Note first that each of the sequences (A n ) and (G n ) consists of IID random elements. 
Let (ri, r 2 , r3) is iV(0, S3) be generic copies of A n and G n . By iterating (13.221 ), we obtain 



@n = Ai(Ai-10n-2 + 0"Gn-l) + oG n 



tn—\ \ n (n—k—1 

Vfc=0 / fc=l \ j=0 



where, as usually, the product is set to be 1 if its index range is empty. We have 

n-l 

n An ~ k = 



k=0 



M n nLi(i+^i,o 



A' 



where 



r n = ea 2 ^A^nW 1 + ^i,fc))- 
j=i fc=i 

Set 5 = max{A, /xE|l + crn|} and note that by (13.131 ) 5 < 1. By triangle inequality and independence of 

ri fc's 



E|T n | < ea 2 ^A n ~ J E 
i=i 



[J(/i(l + ar 1)fc )) 
fc=l 



Similarly, 



^x n E 



ri( i+fTri " 



k=l 



ea 2 ^2\ n - j (/xE|l + ar 1 |) i ~ 1 < e^roS" -1 . 



(/iEll + anl)" < <F' 



It follows that both components of (n^o ^n-kj ©o converge to in probability and thus, this term is 
negligible. 

Since the sequences (A n ) and (G n ) are IID, for every n > 1 we have 

n I n—k—1 \ n I k—1 \ 

e n a^Ag^^IiiaAg,. 

k=l V j=0 / k=l \j=l J 
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By the same argument as above we verify that both components of the sequence of partial sums on the right 
hand side are Cauchy in L\. Hence, the components of the series 



fc-i 



fc=i 



converge in probability (and thus, in distribution). Therefore, the sequence (6 n ) defined by (13.221 ) converges 
in distribution to a random vector X defined in ( I3.24I ). Furthermore, X satisfies the distributional equation 
(13251 □ 



Proof ( Theorem [37Tb : For h = (hi,h 2 ) set Bh := (— oo, hi] x (— oo, h 2 \. Then 

{r = n} = {Qj G B h , j < n, Q n $ B h }, 

so that 

P(r = n) = F(Q n <£ Bhlfy € B h , j < nWQj € B h , j < n) 
= P(^„9 n _i + oG n $ S ft |G n _i G B h )¥(r > n). 

Since O n converge in distribution to X we have 

P(A n G n _! + aG n $ B h , 6 n _i G 



Pn := P(A n G n _i + aG n i B h \Qn-i G 



n— 1 



€B h ) 



:= PjAX + aGtB h XeB h ) as 
^ P(Jf€fl fc ) 

It follows from (13.241 ) that X is symmetric, so since both h\ and h 2 are positive the denominator is at least 
1/2 and does not affect the asymptotics. 

To handle the numerator, using ( 13.261 ), denoting the components of X by X\ and X 2 , and using the 
notation adopted in (13.181) we see that it is equal to 

P((p(l + ar 1 )X 1 + ar 2 , ea 2 X x + XX 2 + ear 3 ) $ B h , (X U X 2 ) G B h ) 
= P(/x(l + ffri)Xi + <rr 2 > /ii, (Xi,X 2 ) G S h ) 
+P(ea 2 Xi + XX 2 + e CT r 3 > h 2) (X U X 2 ) G B fc ) 

-P(/i(l + <7ri)Xi + or 2 > /ii,eo 2 Xi + AX 2 + ecrr 3 > /i 2 ,(Xi,X 2 ) G B h ) 



^^>^^,(X 1 ,X 2 )G^ 
+ p(ra> /t2 " 6a2 f 1 " AX2 ,(X 1 ,X 2 )G^) 



jiXin +r 2 fei - /iXi /» 2 - e«2^i - AX 2 , 
> ,r 3 > , (X 1 ,X 2 )eB h . (3.28) 

Conditionally on (Xl,X 2 ) = (xi,x 2 ), 

/ixiri + r 2 r 3 
Zi := , and z 2 := 
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are N(0, 1) RVs. Hence by letting Fx(x\,X2) denote the distribution function of (Xi, X2), we see that the 
first of the last three probabilities is 

k2 (l - $ ( h-w X) d F x ( Xl ,x 2 ), (3.29) 

Likewise, for the second of these probabilities we get 

We now note that if e is of a smaller order than all other parameters (except possibly a) then (13.101) implies 
that (13.301 ) (and hence also ( 13.281 )) are negligible when compared to ( I3.29I ). To analyze the behavior of ( 13.291 ) 
as a function of its parameters note that by the mean value theorem the quantity in ( 13.291 ) is equal to 



—00 j —00 



for some — 00 < Xq < h. Substituting this into ( 13.271 ) (and neglecting the terms that depend on e) we see 
that 

W(AX + aG<£B h ,XeB h ) f ^ - ^x 

p = — — — — ~ 1 — <£ 



F(XeB h ) V (T<Txa 

If both < /x < 1 and h\ are of order 1 we are in the same situation as with (13- 14b - This shows (13.23I ). □ 



3.5 Diffusive escape 

The exit problems for the stochastic difference equations analyzed in the previous subsections all feature 
the geometric escape mechanism. In the simplest case when the evolution is given by Equation (13 -4b . the 
geometric distribution characterizes the statistics of the times of exit of the trajectories of (13.41 ) from a certain 
neighborhood of the attracting fixed point. In this subsection, we study another important in applications 
statistical regime associated with the exit problem for (13.41 ). the diffusive regime. The role of the diffusive 
regime in characterizing the statistics of the exit times for the trajectories of (13.41 ) is twofold. First, the 
geometric distribution approximates the distribution of the exit times only for sufficiently large times, i.e. 
for large n. In this subsection, we show that in the intermediate range of n, i.e. when n is neither too large 
nor too small, Y n \ are approximated by the sums of the IID RVs and, therefore, the level exceedance times 
are distributed as those for random walks. We refer to this situation as the diffusive regime. Second, we 
recall that to justify the geometric distribution in the proof of Theorem 13.31 we implicitly assumed that the 
rate of attraction of the fixed point is stronger than the noise intensity. Specifically, it is easy to see from the 
proof of Theorem I3.3l that q is required to be o(e), e = 1 — A. The analysis in this subsection does not use 
this assumption. We show that when the noise is stronger than the attraction of the fixed point (albeit both 
are sufficiently small), the mechanism of escape of the trajectories from the basin of attraction of the fixed 
point changes from the geometric to diffusive. Therefore, we conclude this section by pointing out to some 
features intrinsic to the diffusive escape. Specifically, we consider (13.41) . for which as before, we define 

r = inf{/c > 1 : Y k > h}, (3.31) 
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Figure 5: Probability density function corresponding to the distribution ty a (y),a = 1. With a suitable a > 0, 
ty a (y) approximates the distribution of the exit times in the diffusive escape. 

for given h > 0. In contrast to the case considered in Section 3.2, here we assume 

£ = 0(c; a ), a>0. (3.32) 

In Theorem |3.9| below, we show that in the present situation in the intermediate range of n, Y^s behave as 
sums of IID normal RVs. The behavior of the latter is well-known (cf, Lemma |3.11| >. 

Recall that <3?(x) stands for the distribution function of an N(0, 1) RV and denote 

* a (z) = 2 ( 1 - * f J ) , a > 0. (3.33) 



Note that ^ a {x) is a probability distribution function on M + (see Fig. [5]). 



Theorem 3.9. Let the evolution ofY n , n = 0, 1, 2, . . . be given by H3.4\) . Suppose that A = 1 — e with 
e = O (s a ) , a > 0. Then for arbitrary positive (3\ and (32 such that (3\ + 02 < 2a/3, for sufficiently small 
? > 0, 

P( T <n) = ¥ (n)(l + o(l)), a = ^, (3.34) 
in the range ?~' 91 <C n -C ?~ 5 l ~^ 2 . 

Remark 3.10. Since /^i^ > are arbitrary, ^ a (n) practically approximates P(t < n) in the range 1 

n < e~ 2 / 3 . 

We will need the following auxiliary lemma iTTTl Theorem 2.2, Chapter III]. It may be viewed as a 
quantified version of a reflection principle for random walk (see, e.g., l42l Sec. 5.3, 5.4]). 

Lemma 3.11. Let X±,X2, ■ ■ ■ be a sequence of independent, symmetric RVs and set 

k 

Sk = ^^X,-, and S%. = max Sj, j > 1. 

j=i 
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Then for any t,u> the following inequalities hold: 



n 



2F(S n >t + 2u)- 2^F{X k >u)< P(£* > t) < 2F(S n > t). 



(3.35) 



k=l 



Remark 3.12. As was noticed by S. Kwapieri a bit stronger version of the first inequality in (13.351 ) follows 
from a slight modification of the proof of Proposition 1.3.1 in [30]. 

Proof (Theorem I3T91 ) : Without loss of generality, we assume that Yq = (otherwise, apply the same argu- 
ment to Y k — Yq). Note that the distributions of r and Y k * are linked by the following relation 

P(r < n) = F(Y* > h). 

Unwinding (13.41 ) and using Yq = gives 

Y k = ?(A fc_1 n + X k ~ 2 r 2 + ■■■ + Ar fc _! + r k ), 
which we write as S k + W k , where 



We will first show that the main contribution to Y* is from the S 1 ,*. First, by subadditivity of maxima, for 
any < hi < h, 



Further, Y k >S k - \ W k \ so that 

nS* n >h + hi)< F(S* >h+h 1 ,\W n \*<h 1 ) + F(\W n \* > hi) < F(Y* > h) + F(|W„|* > hi), 
which, when combined with (13.371 ) means that 

P(5* >h + hi)- F(\W n \* > hi) < F(Y* >h)< F{S* >h-hi)+ F(\W n \* > hi). (3.38) 
First, we estimate P(| W n \* > hi) in (I3.38I ). To this end, we use 1 - X j = 1 - (1 - e) j < je to obtain 



A- 



fe-1 




(3.36) 



P(^„* >h)< F(S* n + W* > h) < F{S* n >h-hi)+ F(W* > hi) 

< F(S*<h-hi)+F(\W n \*>hi). 



(3.37) 



n-l 




var(W n )= ? 2 ^(l-A^) 2 < ? 2 e 2 



Consequently, by (13.351 ) and (13.36b . we have 



PdWnl* > hi) < 2F(\W n \ > hi) < AF(W n > hi) 
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Next, we turn to estimating the probabilities involving S* in (I3.38I ). By the second inequality in (13.351 ). for 

every u > 0, we have 



F(S* >h-hi)< 2P(S n >h-h 1 )=2P[Z> - jl 1 ) , (3.39) 



while the first one yields 



P(S£>fc + /ii) > 2P (S n >h + hi + 2u)- 2^ P(Vr fc > 



fc=i 



2p fr + fr_+2u\ _ ^ / > _ (3 _ 4Q) 



The combination of (l338i (13391) . and d3~40b yields 

w > ft) > tt (*>*±^)_ w ( z >2)_,r(z> * £), ( , 41) 



P(Y.; >h) < 2? [ Z > !LJ}l) + 4P [ z > A= • ^ | (3.42) 



To complete the proof, we need to chose h\ and u such that 

o(l), -^==o(l), - = o(l), and ^ X ?en 3 / 2 = o(l). (3.43) 



1 + M 

It is straightforward to verify that relations in (13.431 ) hold with hi = ? 2 and u = ? 2 , > 0, 

Pi+/3 2 < 2a/3, and n as in (I3341 . □ 



4 The Poincare map 



In the present section, we consider the type I model, i.e. the randomly perturbed system with the stochastic 
forcing acting via the fast subsystem (see ( 12.91 ) and ( I2.10I )). In the active phase of bursting (when the system 
undergoes spiking), the trajectory of the randomly perturbed system remains in the vicinity of the cylinder 
foliated by the periodic orbits of the fast subsystems, (see Fig. [6^). The time that the trajectory spends near 
L determines the duration of the active phase. The goal of this section is to describe the slow dynamics 
near L. In particular, we will estimate the distribution of the number of spikes in one burst. To this end, 
we introduce a transverse to L crossection S (see Fig. [6^) and construct the first return map. Specifically, 
we estimate the change in the state of the system after one cycle of rotation of the trajectory around L. 
The construction of the first return map for (12.91 ) and (12.101 ) is done in analogy to that for the deterministic 
models of bursting (see |[34l[3TI ). However, the treatment of the randomly perturbed system requires certain 
modifications. First, we have to resolve the ambiguity in the notion of the first return time. The latter 
is due to the fact that generically a trajectory of the randomly perturbed system makes multiple crossings 
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with S during each cycle around L. We refer the reader to the comments following Theorem 2.3 in |fl9l 
for an explicit example illustrating this effect. For the randomly perturbed system, we define the time of 
the first return so that it approaches the first-return time of the underlying deterministic system in the limit 
of vanishing random perturbation. The definition of the first return time motivates the definition of the 
Poincare map (see Definition 4.1). In Sections 4.1 and 4.2, we use asymptotic expansions to construct the 
linear approximation for the Poincare map of the fast subsystem. Here, we use an obvious observation that 
on finite time intervals and for sufficiently small e > 0, the slow variable typically remains in an O(e) 
neighborhood of its initial value. Therefore, for finite times the Poincare map of the fast subsystem captures 
the dynamics of the full system. Since we are interested in long term behavior of the system, to complete 
the description of the first return map we also need to track the (small) changes in the slow variable after 
each cycle of oscillations. This is done in Section 4.3, where we derive a ID map for the slow variable. The 
combination of the ID Poincare map for the fast subsystem and that for the slow variable provides the first 
return map for the full problem ( 12.91 ) and (12- 10b - The linear approximation of the 2D map is used in Section 
4.4 to estimate the distribution of the number of spikes in one burst for the type I model. Effectively, the 
problem is reduced to the exit problem for a ID linear randomly perturbed map. For the latter problem, we 
have already developed necessary analytical tools in Section 3. Finally, in Section 4.5, we comment on the 
straightforward modifications necessary to extend the analysis of this section to cover type II models. 

4.1 Preliminary transformations 

Recall that S stands for the transverse section located as shown schematically in Fig. [6^. Let yo < y^p 

be outside an 0(a) neighborhood of y& p , and xq = (xq,Xq) T G S be from an 0(a) neighborhood of L. 
Consider an initial value problem for (12.91 ) and (12.101 ) with initial data (xq, yo). By standard results from the 
asymptotic theory for randomly perturbed systems |fT9l , we have the following estimate 

yt = y + O(e), (4.1) 

valid on a finite interval of time t G [0, t\. Here and below, for a small parameter e > 0, the symbols O(e) 
and o(e) in the asymptotic expansions of the random functions mean that the corresponding relations hold 
almost surely (a.s.). Specifically, ipt(t) = 0(e) for t £ [tx, t%] means that there exists eo > such that 

sup | e — (e) j < oo a.s.. 

t e [ti,t 2 ] 

e S [0, e ] 

In a similar fashion, we interpret tpt(^) = o(e) when V't(e) is a random function. 
By plugging in (14.11 ) into (12.91 ), we obtain the following SODE 

dx t = f (x t ) dt + ap(x t )dw t + O(e), (4.2) 

where / (x) := / (x, yo), p(x) := p(x,yo) , and yo is fixed. Equation (14.21 ) with e = a = has an 
exponentially orbitally stable periodic solution x = <ft(t, yo) of period T(yo): 

L(y ) = {x = </,(9, yo): G [0, T(y ))} (cf. Q). 
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To simplify the notation, throughout the analysis of the fast subsystem, we will omit to indicate the depen- 
dence on uq when refer to L, (j), and T. At each point x = <p(9) G L, we define vectors 

r(9) = (f\x),f 2 (x)) T and u{9) = Jf(x), where J= ( J ~ 1 V (4.3) 

pointing in the tangential and normal directions, respectively. To study the trajectories of (14.21 ) in a small 
neighborhood of L, it is convenient to rewrite (14.21) in normal coordinates (9, £) [23 ]: 

* = + 0€[O,T). (4.4) 



Lemma 4.1. For sufficiently small 5 > Equation ( 14.41 ) defines a smooth change of coordinates in 

B 5 = {x = cj>(9) + fr{9) : |e| < 5, G [0, T)}. 

/n new coordinates, ( 14. 2D /zo^ the following form: 

M t = (l + b 1 (9 t )^)dt + ah 1 (9 t ^ t )(l + b 2 (9 t )Ct)dw t + 0(e,5 2 ,a 2 ), 
<%t = a(9 t )t t dt + o-h 2 (9 t ,t t )dw t + 0(e,5 2 ,a 2 ), 

where smooth functions a(9), b\(9), and b 2 {9) are T —periodic and 

0< // := exp ^jf a(0)d0^ = exp f / div/ (0(0))^) < 1, 
1 f l , „2 f 2 







= = — r^2 — ' w>£) 



< r,r > 



i/r 



< r,r > 



i/r 



(4.5) 



(4.6) 
(4.7) 



(4.8) 
(4.9) 



Proof : The proof of the lemma follows the lines of the proof of Theorem VI. 1.2 in E31 . Let z 
(z 1 , z 2 ) T := (9, £) T and denote the transformation in ( 14.41 ) by 



Note 



\Dv( 



x = v(z), z G Bs- 

|/(^))|Vo, 0g [0,T). 



(4.10) 



<p\0) -f 2 {4>{9)) 
2 '(9) f 1 (ct>(9)) 

Therefore, for sufficiently small 8 > 0, (14.101 ) defines a smooth invertible transformation in Bs- Denote the 
inverse of v by z = u(x), x G v (Bs) and note that 



By Ito's formula, we have 
and, therefore, 



[Dn(x)]- 1 = Dv(z), x G v(B s ). 

dz t = Du(xt)dxt + 0(a 2 )dt 
Dv(z t )dz t = dx t + 0(a 2 )dt. 



(4.11) 

(4.12) 
(4.13) 
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Figure 6: (a) Crossection S is used in the construction of the first return map. (b) The phase plane of the 
fast subsystem d2.il ) for y G (y S n, Vbp)- 



24 



By recalling that z = (6, £) and after plugging in (14.21) into (14.131) . we obtain 

dB t + v(e t )d£ t = UWt)) + Df{<l>(0 t ))v(0t% + Q{0tAt))dt 



H 777- £t 
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+ apdw t + 0(e, cr 2 ), 



(4.14) 



where 



Note that 



,o = f one) + - / w)) - Df (m) *>m = ° (£ 2 ) , iei < S. 



f (m) = t(o), t t (9)t(0) = v(e) T v{o) = i/ (mtf , 
d -jf(m) = JDf(m)f(m)- 



dv{9) 



(4.15) 
(4.16) 



d9 d9 

Taking into account (14.151 ) and (14.16I ). we project (14.141 ) onto the subspace spanned by T(9 t ) and after some 
algebra obtain: 

f T Q + f [DfJ - JDf] fi t + af T pw t + 0(e) 



1 ' ' Ff + FJDfftt 

Here and for the rest of the proof, for brevity we use the following notation: 

f := f (<!>&)) , Q:=Q(e t ,£ t ), and u := v{6 t ) 
Equation ( 14.171 ) can be rewritten as (14.61 ) with 

1 

T7P 

b 2 (9t) = 



(4.17) 



h(9 t ) 



,f T [DfJ -JDf] f, 



i 



I/I' 



f T JDff. 



Similarly, by projecting (14.141 ) onto the subspace spanned by v(9) and using (14.151 ) and (14.141 ). we derive 

it = a(6 t )t t + oh 2 (6 t )w t + 0(5 2 ) 1 

where 



a(9 t 



1 



Dfu + 



dv 
d9 



2v T dv 



[v[ 2 d9 



The expression in the square brackets can be simplified as follows: 

dv 



Dfv + 



d,e 



f [J T DfJ + Df] =div/(<X0))|/| 2 . 



Also, 



Therefore, 



2v T dv 2 T T d 

tJ J d~9 J 



\v\ 2 d9 \f\< 



\ff d9 J \f\*dB 



l -A\n 2 = 4- Mf m))\ 2 ■ 



d9 



a(0) = dhr/ (*(*)) -^ In |/(W))| 2 . 



(4.18) 



Equation (14.181 ) implies (14.81) . since the integral over [0, T] of the last term on the right hand side of (14.181 ) 
is zero. □ 
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4.2 The Poincare map for the fast subsystem 



In the present subsection, we analyze the trajectories of the randomly perturbed system (14.21 ) lying close to 
the limit cycle L(yo)> Uo < Vbp- To this end, we consider an IVP for (14.61 ) and (14.71) subject to the initial 
condition: 

9 = and |£ | < S. (4.19) 

Throughout this section, we assume (even when it is not stated explicitly) that 5 > is sufficiently small. It 
will be convenient to view the range of 6 t as R 1 rather than a circle. Equation ( 14.41 ) provides the transforma- 
tion of {9 t , £t) to the Cartesian coordinates even when 6 t exceeds T. 

We now turn to the construction of the Poincare map. Condition 6 = defines a transverse crossection 
of L(yo), S. The trajectory of the deterministic system (14.61 ) and ( 14.71 ) with a = returns to S in time 
T + O(£o)- To define the Poincare map for the randomly perturbed system, we also use another transverse 
crossection E, which is located at an 0(1) distance away from S. Let (9t,£,t) be the solution of the IVP 
(|43T> . (14771) . and (I4T91 and 

T = inf{i>0: (0t,&)e£}. 

Definition 4.2. By the time of the first return of the trajectory ( 14.61 ), ( 14. 71 ), and < \4.19i to "£, we call stopping 
time T such that 

T = mf{t >f: 9 t = T}. (4.20) 
The first return map for (14. 61 ), (14. 71 ), and A4.19i is defined as 

£ = P(€o), where £ = £t- 



In the remainder of this subsection, we compute the linear part of the Poincare map. In the asymptotic 
expansions below, we omit to indicate the dependence of the remainder terms on e > 0. The latter is 
assumed to be sufficiently small so that it has no effect on the leading order approximation of the Poincare 
map. 

The following notation is reserved for four functions, which will appear frequently in the asymptotic 
expansions below: 

A(t, s) = exp{ J s * a(u)du}, A(t) = A(t, 0), 
B(t,s) = fl A(u, s)bi(u)du, B(t) = B(t,0). 

Lemma 4.3. On a finite time interval t G [0,t], < i < oo, the solution of the IVP ( 14.61 ), ( 14. 71 ) and ( 14. 1 91 ) 
admits the following asymptotic expansion 

6 t = e^+aO^ + Oia 2 ,^), (4-21) 
it = d 0) + < (1) +O(a 2 ,£ 2 ). (4.22) 

The leading order coefficients are given by 

ef ] = t + ^B(t) + O(f ), (4.23) 
e t (0) = £oA(t)+0(e o ). (4.24) 
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The first order terms are given by Gaussian diffusion process Zt = ( 6^ , ^ 



zt 

>0 



fu(t, s)h(s)dw s + 0(Co), (4.25) 
Jo 



where 



n/ -- s) 5 A(t,s) ) ' h(t) := Ht > 0) = ( h i(t> Q )Mt,0)f • ( 4 - 26 ) 



Proof: The procedure for constructing asymptotic expansions of solutions for a class of IVP, which includes 
(14.61) . (14.71) and ( 14.191 ) can be found in (HQS!. These sources also contain the estimates controlling the 



remainder terms. The coefficients df* and £,1°^ are determined as follows. By plugging in (14.211) and 
(14.221 ) into (14.61 ) and (14.71) and extracting the coefficients multiplying different powers of a, one obtains 
IVPs for the functions on the right hand sides of (14.211 ) and (14.22b - Specifically, for the leading order terms 
we have the following IVP: 

0(°) _ i , fa(0)\ c(°) 



i+6x cur. ( 4 - 27 ) 



To the next order, 



4 (0) = a(e^)tf\ (4.28) 
4 0) = Co, ef ] = 0. (4.29) 



A(t,£o)z t + h (0 t (O) ,£ t (O) W, (4.30) 



= 0, (4.31) 



where = (6> t (0) , ) , ft = (hi, h 2 ) T , and 




Here, we explicitly indicated the dependence of the leading order coefficients on £o an d used prime to 
denote the differentiation with respect to 9. Formulae (14.23l) - (14.26l) in the statement of the lemma follow 
from ( |4.27b -( l4T32b . The details can be found in the appendix to this paper. 

□ 

Next, we calculate the time of the first return. 
Lemma 4.4. The time of the first return is given by 

T = T (0) + aT (1) + o(a) + O(£o), ( 4 -33) 

where 

T<® = T-Z Q B{T) + 0{Zl), (4.34) 

T (1) = -a9>p = -o- [h 1 (u) + B(T,u)h 2 (u)]dw u . (4.35) 

Jo 
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Proof: From the definition of the first return time, (14.211) . and (14.23b . we have 

T + CoB{T) + aO^ + 0{o 2 ^l) = T a.s.. (4.36) 

Thus, 

lim T = T(°)(£o) a.s., (4.37) 

cr^O 

where T(°'(£ ) is found from the following equation 

T (0) (&) + (V (0) (Co)) + O(e 2 ) = T. (4.38) 



Equation ( 14.381 ) implies (14.341 ). Furthermore, the combination of (14.341 ), ( 14.361 ), and (14.371 ) yields (I4.35D . 

□ 



Lemma 4.5. The first return map is given by the 

£ = (1 + + ^ 2 + ((t) + O(£ 2 ), (4.39) 
where Gaussian RVs r\ >2 are given by 

n = -a(0) I [hi(u) + B(T,u)h 2 (u)]dw u , r 2 = I A(T,u)h 2 (u)dw u . (4.40) 
Jo Jo 

Proof: From (l4~22l , ( |4T24T )-( |4~26T ), and (14331 , we have 

f = = toA(T) + a [ T A(T, s)h 2 (s)dw s + 0(a 2 

Jo 

= ^A(T)+ar 2 + O(a 2 ,e ), (4.41) 
where r 2 is defined in (14.40b . The first term on the right hand side of (14.411 ) can be rewritten as follows 

A(T) = A{T)A(T + aT ( - 1 \T) + o(a) + O(^) = ^exp(aa(0)T ( - 1 A+o(a) 



= fi(l- aa{0)e ( T ] ) + o(a) + O(fo). (4.42) 
Finally, we extract the expression for o!p from ( 14.251 ) and ( 14.261 ): 

9^ = / [h^u) + B(T,u)h 2 {u)]dw u . (4.43) 
Jo 

Equations (@gl) - (@g3) yield (14391 and (14401 . □ 



Remark 4.6. We close this section by observing that as follows from (14.401 ) RV r\ and r 2 are stochastic 
integrals of different deterministic functions, say f(t) and g(t) with respect to the same Brownian motion 
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over the interval [0, T]. Consequently, their joint distribution is bivariate normal with mean vector and a 
covariance matrix that whose diagonal entries are 

f 2 (t)dt and f g 2 (t)dt, 
o Jo 



and the off diagonal entry is 



f(t)g(t)dt. 



o 



This is perhaps easiest to see by using Riemann representation of a stochastic integral (see e.g. B2l Propo- 
sition 7.6]), basic properties of Brownian motion, and a fact that a random vector is multivariate normal if 
and only if any linear combination of its components is a normal RV. 



4.3 The first return map for the slow variable 

Our next goal is to estimate the change of the slow variable, y t , after one cycle of oscillations of the fast 
subsystem for the following initial conditions: 

< y bp ~ yo = 0(1), x = 0(0) + £oKO) G E, and |£ | < 6. (4.44) 

We denote the first return map for y by 

V = P(y, Co), where P(y , f ) = y T , 

and T is the first return time of the fast subsystem (see Defmition |4.2| ). 

Lemma 4.7. The first return map for y has the following form: 

P(y, = V + tG{y) + ear 3 + eat + o(ea), (4.45) 

where ^ 

G(y) = [ g^(s),y)ds (4.46) 
Jo 

and r% = N (0, 0(1)) and a is a constant independent of a and e. 

Remark 4.8. Recall that T and (/>(•) are functions of slow variable y (see (12.21) ). To avoid using cumbersome 
notation we continue to suppress the dependence on y. 

Proof: Bv dlTOl . 

fT 

yr = yo + e g(x s , y )ds + 0(e 2 ), (4.47) 
Jo 

where x s satisfies IVP (l4~6l) . (l4~7i and (14- 19b . Let x = (f)(9) + £i/(0) and denote 

g(9,Z,y):=g(x,y), g (s) = ~g(s, 0), 9l (s) = ||( S) 0), and 52 ( 5 ) = ||( 5j 0). (4.48) 
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Using (14.481 ). we rewrite (14.471 ) as 

vt = y + e /%(^ 0) + ^ 1} ,d 0) + + 0(6a 2 

JO 



Using the Taylor expansion for g in ( 14.491 ) and (14-2 lb . (14.221) and (14.331) . from ( 14.491 ) we derive 

rT 



VT 



yo 



+ 



'0 



+ e J {g (s) + 9l (s) [toB(s) + a0«] + g 2 (s) [^ A(s) + } 

gd(a)ds + o(ea) + 0(efg). 
We approximate the last integral on the right hand side of (14.501 ) by 

g ( S )ds = - 50 (0) k B(T) + <rfl«l + o(<r,£ ). 



(is 



r 



T 



The combination of (14.501 ) and ( |4.5 lb implies d4.45b with 



fcW^W + 92(s)A(s)} ds - g (0)B(T), 



T3 



gi(s)6W+g 2 (s)^ 



ds. 



(4.49) 



(4.50) 



(4.51) 



(4.52) 

(4.53) 
□ 



4.4 The exit problem 

In the present subsection, we first combine the return maps derived for the slow and fast subsystems to 
obtain the Poincare map for the full three-dimensional system. Next, we approximate the Poincare map 
and the BA of the limit cycle L(y c ) and characterize the distribution of the exit times for the approximate 
problem. This distribution is then related to the distribution of the number of spikes within bursting episodes. 
To approximate the Poincare map we linearize it around the stable fixed point of the deterministic map 
corresponding to the limit cycle L(y c ). Aside from the systematic derivation of the Poincare map in the 
previous subsections, we offer no rigorous justification for substituting the nonlinear Poincare map with its 
linear part in the analysis of the exit problem. While in general, such approximation may not be accurate, 
we believe that for the present problem, the analysis of the linearized system captures the statistics of the 
first exit times well for the following reason. In models of square wave bursting the limit cycle generating 
spiking is often located close to the boundary of its BA (see Fig. [6}) for a representative example). Therefore, 
before the trajectories leave the BA, they remain in a small neighborhood of the limit cycle, where the linear 
part of the vector field governs the dynamics. After these preliminary remarks, we turn to the derivation of 
the approximate problem and its analysis. 
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Lemmas 14.51 and 14.71 yield the asymptotic formulae for the first return map of the randomly perturbed 
system ( 12.91 ) and (12.101) in the normal coordinates (14.41) : 



£ n+1 = iiin (1 + ar 1>n ) + ar 2 , n + o(a), (4.54) 
Vn+i = Vn + eG(y„) + eoT 3) „ + ea£ n + o(ecr), n = 0, 1, 2, . . . , (4.55) 

where (£q, yo) are given in (14.441 ) and the expressions for a and rj jTl , i = 1, 2, 3 are are given in (I4.40b . (l4.52b . 
and (I4.53I ). Recall that by (SS) (see Section 2), G(y) has a simple zero at y = y c and A := —G'(y c ) > 
0. Thus, (0, y c ) is an attracting fixed point of the unperturbed map ( 14.541 ) and ( 14.551 ) with a = 0. The 
linearization of ( 14.541 ) and ( 14.551 ) about (0, y c ) yields 

£ n+ i = fj£ n (1 + ah,n) + crr 2 , n , (4.56) 
Vn+l = Ar?„ + eaf 3tn + ea 2 £„, ra = 0,1,2,..., (4.57) 

where f] = y — y c , 0<A = 1 — eai, and < \i < 1. The distributions of the RVs r$ n , i = 1,2,3 
depend on y n , as both the upper bound of integration T and the integrands in ( 14.401 ) and ( 14.531 ) are smooth 
functions of y. The stochastic terms r\ n , i = 1,2,3 in the linearized system are obtained by evaluating 
the expressions for fj )Tl , i = 1, 2, 3 in (14.401) and (14.531 ) at y = y c . Thus, (fi . n , ^2,n> ^3,n) are HD copies of 
a (0, S3), where the entries of S 3 are O(l) in a sense that they do not depend on any other parameters. 
Further, we approximate the BA of L(y c ) by a cylindrical shell, so that in (£, 77) coordinate plane, it projects 

to n := —h£,h,£ x —hrjjhrj for some h^„ > > independent of cr > 0. Each iteration of the 
Poincare map corresponds to a spike within a burst. The burst terminates when the trajectory leaves the BA 
of L(y c ). Assuming that the linearization (14.561 ) and (14.571) and II provide suitable approximations for the 
Poincare map and the BA of L(y c ) respectively, the distribution of the number of spikes in one burst can be 
approximated by the distribution of the first exit times for the trajectories of (14.561 ) and (14.571) from II: 

t = min{T£, r v }, (4.58) 

where 

-it = inf {£ n > he} and r„ = inf {q n > h n }. 

n>0 n>0 

We are now in a position to apply the the results of Section 3 to describe the distribution of (I4.58I ). By 
Theorem l3.7l the distribution of r is asymptotically geometric with parameter 

a c 

(4.59) 



CV2TT 

for some C > independent of e and a. In the proof of Theorem 13.71 we studied a class of 2D randomly 
perturbed maps that includes (14.561) and (14.57b - However, the distribution of r is effectively determined by 
the first equation (14.56b . i.e. by the ID first return map of the fast subsystem. This can be seen by observing 
that according to the approximations given at the end of proof of Theorem l3.7l (see the arguments following 
(13.301 )) if e > is sufficiently small then tc -C t v and r ~ 7£. Thus, in type I models the distribution of 
spikes in one burst is effectively determined by the ID first return map for the fast subsystem (14.56b . In 
particular, the statistics of the number of spikes in one burst does not depend on the relaxation parameter 
e > 0, provided the latter is sufficiently small. 
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4.5 Type II model 



The derivation of the Poincare map for the type II models differs from the analysis in Sections 4.1-4.4 for 
type I models only in some minor details. In this subsection, we comment on the necessary modifications 
and state the final result. Recall that in contrast to type I models, in ( 12.111 ) and (12.121 ), stochastic forcing 
enters the slow equation. As before, the initial condition is given by (14.44b - On finite time intervals, solutions 
of the IVP for (12.111) and (12.121) admit the following asymptotic expansions 

x t = xf ] + eaxf ] + O ((ea) 2 ) , (4.60) 
Vt = yl° ] + eay\ l) + O ((ea) 2 ) . (4.61) 

where the first order corrections xf^ and y^ are Gaussian processes (cf. Theorem 2.2 Ifl9l0 . Using (14.601 ) 
and (I4.61I ). we obtain the leading order approximation of the fast subsystem: 

ff \ , d f( x ?\yo) (l) , / n fA ™x 

xt = f(xt,Vo) + ecr — y\ + o{ea). (4.62) 

From this point, the derivation of the Poincare map follows the same lines as we described in detail for type 
I models in Sections 4.1-4.4. We omit any further details and state the final result, the linear approximation 
of the Poincare map for the present case: 

£n+l = fJ-Cn (1 + £(7h,n) + ^h,n, (4.63) 
r] n+1 = Xr] n + eaf^ n + ea 2 £n, n = 0,1,2,..., (4.64) 

As in the previous case, we are interested in the distribution of the first exit time r (see ( I4.58I )). To estimate 
the latter, we use the same argument as in the previous subsection. This time the system is described by 

6 n+ i = A n+1 @ n + creG n+ i, n > 1, (4.65) 
. The presence of the factor e in both components of G n leads to 



where A„ is as before and G r , 



r2,n 

the following expression for the numerator of p (see (13.271 )): 



(MXirx + n hi-fiXx h 2 -a 2 Xi X(h 2 - X 2 ) , v v ,^ v 
> ,^3 > 1 , (Ai, X 2 ) G k>h 

This expression decays very fast as a function of h 2 — X 2 and since X 2 has heavy tails it is approximated 
(up to inessential polynomial factors) by 



aXxn + r 2 hi - /iXi h 2 
> ,r 3 > 



0Xx to-axx cr 



^,(x 1 ,x 2 )eB h y 



We are now in the analogous situation to that encountered in (I3.28I ). except that the small parameter e > 
appears in the denominator of the other variable. As a consequence, this time we obtain that <C t v for 
small e > 0. Therefore, in contrast to type I models, the escape of a trajectory of (12.111 ) and (12.121 ) from A 
is dominated by the slow subsystem, i.e., r = t v . 
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5 Numerical example 



In the present section, we illustrate the statistical regimes identified in this study with numerical simulations 
of a conductance based model of a neuron in the presence of noise. To this end, we use a three variable 
model of a bursting neuron introduced by Izhikevich f26l| . The model dynamics is driven by the interplay 
of the three ionic currents: persistent sodium, Ino,p, the delayed rectifier, Ik, a slow potassium M-current, 
1km, and a passive leak current 1^. The following system of three differential equations describes the 
dynamics of the membrane potential, v, and two gating variables n and y: 

Cv = F(v,n,y), (5.1) 
T n n = raoo(w) - n, (5.2) 
T y y = yoo{v) - y, (5.3) 

where F(v, n, y) = -gNapm^v) (v - E Na p) - g K n(v -E K )~ 9kmv{v - E K ) -g L (v-E L )+I; g s and 
E s , are the maximal conductance and the reversal potential of I s , s G {NaP, K, KM, L], respectively; and 
/ is the applied current. The time constants r n and r y determine the rates of activation in the populations 

of K and KM channels. The steady-state functions are defined by s^v) = fl + exp (tt^)) > s e 
{m, n, y} . The parameter values are given in the caption to Fig. |7] This completes the description of the 
deterministic model. The random perturbation is used in the form of white noise, awt and is added to the 
first equation (15.11 ) for type I model or to the third one (15.31 ) for type II model. After suitable rescaling, 
these models can be put in the nondimensional form ( 12.91 ), (12.101) or (12.1 II) . (12.121) . The separation of the 
timescales in the nondimesional models (i.e. small e >)) is the result of the presence of the disparate time 
constants 17, 3> r n in the original model (see caption to Fig. |7J). 

The parameters of the deterministic system are chosen so that it has a limit cycle located as shown in 
Fig. [3fc. In the presence of small noise the system generates bursting. In each numerical experiment, we 
integrated the randomly perturbed system using the Euler-Maruyama method |[24ll until it generated 5, 000 
bursts. We used these data to estimate the probability density for the number of spikes within one burst. 
In Fig. 13 we plot the histograms for the number of spikes in one burst for type I and type II models. 
The histograms in Fig. [7] are scaled to approximate the probability density function (PDF) for the number 
of spikes in one burst. Both PDFs shown Fig. |7^,b have distinct exponential tails as expected for the 
asymptotically geometric RVs. Note that the distribution in Fig. fits well with the geometric distribution 
for N > 100, while in Fig. the geometric distribution fits the data almost on the entire domain N > 10. 
In addition, the peak in the histogram in Fig. is reminiscent of the PDF characteristic for the diffusive 
escape (see Fig.[5]>. For comparison, we plotted a slightly shifted diffusive PDF, ty a (x), a = 10.8 in Fig.|7^. 
Matching the data and ^ a is a delicate matter, because it is not clear how wide is the range of n, to which 
the estimates of Theorem 13.91 apply. Nonetheless, the qualitative similarity of the peak in the histogram in 
the range n ~ 50 — 100 and the diffusive PDF is apparent. We repeated these numerical experiments for a 
few other sets of parameters and found qualitatively similar results. 

Collecting the statistical data shown in Fig. [7] requires integrating the system over very long intervals 
of time, for which it would be hard to justify the accuracy of the Euler-Maruyama method. However, 
capturing the statistical features of the dynamical patterns does not require having an accurate solution on 
the entire interval of time, because they are determined by the discrete dynamics of the first return map. The 
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Figure 7: The histograms for the number of spikes in one burst. The histograms computed for the type I 
model in (a) and type II in (b) are normalized to approximate the corresponding PDFs. The tails of both 
functions are well approximated by the exponential densities with parameters 0.0067 and 0.0125 respec- 
tively. In (b) the exponential distribution already gives a very good approximation for the number of spikes 
exceeding 10. The region of exponential behavior in (a) starts around n ~ 100. In (a), we also plotted in 
solid blue line the shifted diffusive density ^ a {x — 25), a 10.8. Although it is hard to claim a quantitative 
fit of the diffusive density and the data, the qualitative similarity between the diffusive pdf ^> a (x) and the 
peak in the data in the range n ~ 50 — 100 is apparent. The values of parameters are C = 1 (^Fcm~ 2 ) ; 
9na = 20, g K = 10, g KM = 5, g L = 8 (mScm' 2 ); E Na = 60, E K = -90, E L = -80 (mV); 
a m = —20, a n = —25, a y = —10 (mV); b m = 15, b n = 5, b y = 5; r n = 0.152, r y = 20 (ms^ 1 ), 
/ = 5pA, and a = 1. 
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iterations of the latter are expected to be insensitive to the numerical noise as suggested by the analysis of 
the randomly perturbed maps in Section 3. Therefore, we only need to have accurate numerical solutions on 
the time intervals comparable with the typical periods of the fast oscillations. This is easy to achieve with 
the Euler-Maruyama method. We repeated these numerical experiments using the second order Runge-Kutta 
method and obtained very similar results. These informal arguments form the rationale for using the above 
numerical scheme. The rigorous justification of the numerics is beyond the scope of this paper. 

Acknowledgments. This work was partially supported by NSF grant IOB 0417624 (to GM) and NSA grant 
MSPF-04G-054 (to PH). 



Apendix 

In this appendix, we provide the details of the derivation of (14.23 l )- (14.26l ). which were omitted in the main 
part of the paper. 

To derive (14.231 ) and (14.241 ). we first note that Q±' is a monotonic function on [0,t\, provided 5 > is 
sufficiently small. Thus, 



a 



(9^+0(3), 



d0(°) 
and 

^)(0(O) ) = ^ A{ 0(O) ) + O (f o ). (A.1) 
By plugging in (1A.1I) into (14.271 ). we have 

0(°) = l + 6l (0(°) )£ A(P<®). (A.2) 

By Gronwall's inequality, 

t (o) =^ + O(£g), *e [0,t\, (A.3) 

where ip t solves 

4 0) =l + Z b l (t)A(t), Vo = 0. (A.4) 
The combination of (1A.1I) . (1A.3I) . and (1A.4I) implies (14.241) . 

We next turn to IVP (14.301) . (14.311 ) and (14.241) . Let U(t, £o) denote the principal matrix solution of the 
homogeneous system 

zt = A(t,£ )z t . (A.5) 
Then the solution of ( 14.301 ) and ( 14.311 ) is given by 

zt = J* U(t, s, £ )h (9f\ dw s = J U(t, s)h (s, 0) dw s + O(£o), t G [0, t\, (A.6) 

where 

U(t,s,^) = U(t,^)U- 1 (s,^) and U(t, s) = U(t, s, 0). (A.7) 
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Finally, by integrating (1A.5I) with £o = and appropriate initial conditions, one computes 

*m>- (;$>). 

The expression for U(t, s) in (14.261 ) follows from (IA.7I ) and (1A.8I) . 
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